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Abstract One main issue, when numerically integrating autonomous Hamiltonian 
systems, is the long-term conservation of some of its invariants, among which the 
Hamiltonian function itself. For example, it is well known that classical symplectic 
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a new family of methods, called Hamiltonian Boundary Value Methods (HBVMs), 
is introduced and analyzed. HBVMs are able to exactly preserve, in the discrete so- 
lution, Hamiltonian functions of polynomial type of arbitrarily high degree. These 
methods turn out to be symmetric, precisely A-stable, and can have arbitrarily high 
order. A few numerical tests confirm the theoretical results. 
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The numerical solution of Hamiltonian problems is a relevant issue of investiga- 
tion since many years: we refer to the recent monographs ll8l [T4l for a comprehensive 
description of this topic, and to the references therein. 

In a certain sense, the use of a numerical method acts as introducing a small per- 
turbation in the original system which, in general, destroys all of its first integrals. 
The study of the preservation of invariant tori in the phase space of nearly integrable 
Hamiltonian systems has been a central theme in the research since the pioneering 
work of Poincare, the final goal being to asses the stability of the solar system. From 
a numerical point of view, results in this respect are still poor, and this is justified 
by considering the delicacy of the problem: as testified by KAM theory, even small 
Hamiltonian perturbations of completely integrable systems, do not prevent the dis- 
appearance of most of the tori, unless a Diophantine condition on the frequencies of 
the unperturbed system is satisfied. 

At the times when research on this topic was started, there were no available 
numerical methods possessing such conservation features. A main approach to the 
problem was the devising of symplectic methods. However, though the numerical so- 
lution generated by symplectic (and/or symmetric) methods shows some interesting 
long-time behavior (see, for example, (8j Theorems X.3.1 and XI. 3.1]), it was ob- 
served that symplecticity alone can only assure, at most, the conservation of quadratic 
Hamiltonian functions, unless they are coupled with some projection procedure. In 
the general case, conservation cannot be assured, even though a quasi-preservation 
can be expected for reversible problems, when symmetric methods are used (see, 
e.g., ID). On the other hand, a numerical "drift" can be sometimes observed in the 
discrete solution 0. One of the first successful attempts to solve the problem of loss 
of conservation of the Hamiltonian function by the numerical solution, is represented 
by discrete gradient methods (see 031 and references therein). Purely algebraic ap- 
proaches have been also introduced (see, e.g., J6)), without presenting any energy- 
preserving method. 

A further approach was considered in fl6l , where the averaged vector field method 
was proposed and shown to conserve the energy function of canonical Hamiltonian 
systems. As was recently outlined (see Q), approximating the integral appearing 
in such method by means of a quadrature formula (based upon polynomial inter- 
polation) yields a family of second order Runge-Kutta methods. These latter meth- 
ods represent an instance of energy-preserving Runge-Kutta methods for polynomial 
Hamiltonian problems: their first appearance may be found in ifTol . under the name of 
s-stage trapezoidal methods. Additional examples of fourth and sixth-order Runge- 
Kutta methods were presented in ifTTl and ifTJl . 

In H10II1 ll[T3l . the derivation of such energy preserving Runge-Kutta formulae 
relies on the definition of the so called "discrete line integral", first introduced in 
m~2l . However, a comprehensive analysis of such methods has not been carried out so 
far, so that their properties were not known and, moreover, their practical construction 
was difficult. 

In this paper we provide such an analysis, which allows us to derive symmetric 
methods, of arbitrarily high order, able to preserve Hamiltonian functions of poly- 
nomial type, of any specified degree. Such methods are here named Hamiltonian 
Boundary Value Methods (HBVMs), since the above approach has been, at first, stud- 
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ied (see, e.g., 11 1U13I0 in the framework of block Boundary Value Methods. The latter 
are block one-step methods (5). However, for the sake of clarity, and later reference, 
the equivalent Runge-Kutta formulation of HBVMs will be here also considered. 

In the remaining part of this section, we introduce the background information 
concerning the approach. Let then 

y'=JVH(y), y(0)=y eR 2m , (1.1) 

be a Hamiltonian problem in canonical form, where, by setting as usual /,„ the identity 
matrix of dimension m, 




(1.2) 



and where the Hamiltonian function, H{y), is a polynomial of degree v. It is well 
known that, for any y* e K 2m , 

H(y*)-H(y )= [ VH(y) T dy = [ a'(t) T VH(a(t))dt, (1.3) 

Jyo-tf JO 

where a : [0,1] — ► R 2m is any smooth function such that 

ff(0)=yo, cr(l)=y*. 
In particular, over a trajectory, y(t), of (II. lb . one has 

H(y(t))-H(yo) = f VH(y(x)) T y' {x)d% 
Jo 

= fvH{y{T)) T JVH{y(T))dT = Q, 
Jo 

due to the fact that matrix J in ( 11.21 ) is skew-symmetric. 

Here we consider the case where (7(f) is a polynomial of degree s yielding an 
approximation to the true solution y(t) in the time interval [0,h] which, without loss 
of generality, is hereafter normalized to [0, 1]. More specifically, given the s+ 1 ab- 
scissae 

= co<ci < ■■■ <c s = 1, (1.4) 

and the approximations y; w y(ci), is meant to be defined by the interpolation 
conditions 

a(ci)=y u ; = 0,...,s. (1.5) 

Actually, the approximations {y,} will be unknown, until the new methods will be 
fully derived. 

A different, though related concept, is that of collocating polynomial for the prob- 
lem, at the abscissae (11.41 ). Such a polynomial is the unique polynomial u(t), of degree 
s+ 1, satisfying 

u(co)=yo, and «'(c,-) = JVH(u(cj)), i = 0,...,s. (1.6) 
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It is well known that (11.6b define a Runge-Kutta collocation method. Moreover, 
the set of abscissae ( 11.41 defines a corresponding quadrature formula with weights 



h '= II — —. dt > i = 0,l,...,s, (1.7) 



./ <>../''<'' 

which has degree of precision ranging from s to 2s — 1, depending on the choice of 
the abscissae (11 .41 . In particular, the highest precision degree is obtained by using the 
Lobatto abscissae, which we shall consider in the sequel. [j The underlying collocation 
method has, then, order 2s. 

Remark 1 Choosing a Gauss distribution of the abscissae {c, } raises the degree of 
precision of the related quadrature formula to 2s + 1. In such a case, it is interesting 
to observe that applying ( 17.3P along the trajectory u(t) and exploiting the collocation 
conditions il.6V one gets 

H(u(c s ))-H(y ) = [ u'(t) T VH(u(t))dt (1.8) 
Jo 

s 

= Y,b i u'(c i ) T VH(u(c i ))+R s = R s , 

i=0 

where R s is the error in the approximation of the line integral. Therefore, H(u(c s )) — 
H(yo) if and only if R s = 0, which is implied by assuming that the quadrature for- 
mula with abscissae {c,} and weights {bf\ is exact when applied to the integrand 
u'(t) T VH(u(t)). However, since the integrand has degree 

s+(v-l)(s+l) = v(s+l)-l, 

it follows that the maximum allowed value for V is 2. Indeed, it is well known that 
quadratic invariants are preserved by symmetric collocation methods. On the other 
hand, when V > 2, in general R s does not vanish, so that H(u(c s )) ^ H(yo). 

The above remark gives us a hint about how to approach the problem. Note that 
in ( 11.8b demanding that each term of the sum representing the quadrature formula is 
null (i.e., the conditions dl.6l )). is an excessive requirement to obtain the conservation 
property, which causes the observed low degree of precision. A weaker assumption, 
that would leave the result unchanged, is to relax conditions (11.61 ) so as to devise a 
method whose induced quadrature formula, evaluated on a suitable line integral that 
links two successive points of the numerical solution, is exact and, at the same time, 
makes the corresponding sum vanish, without requiring that each term is zero0 

If we use <j(t) instead of u(t ), the integrand function in (11.31 ) has degree vs — 1 
so that, in order for the quadrature formula to be exact, one would need say, k + 1 
points, where 

t1- 



Different choices of the abscissae will be the subject of future investigations. 
2 More precisely, in the new methods, conditions U.6t will turn out to be replaced by relations of the 
form o"'(c,) = Y,j fiijJ^H(o(cj)), which resemble a sort of extended collocation condition (see also [13] 
Section 2]) since o"'(c,) brings information from the global behavior of the problem in the time interval 
[0, ft] (see )3.U - f3~8l in Section|3]and the analogues in Section[4}- 
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if the corresponding Lobatto abscissae are used. Of course, in such a case, the vanish- 
ing of the quadrature formula is no longer guaranteed by conditions ( 11.6b and must 
be obtained by a different approach. For this purpose, let 

r = k-s, (1.10) 

be the number of the required additional points, and let 

0< Ti < ••• < X r < 1, (1.11) 

be r additional abscissae distinct from dl.4l ). Moreover, let us define the following 
silent stages |fl~3), 

Wi = a{Xi), i=\,...,r. (1.12) 

Consequently, the polynomial o(t), which interpolates the couples (c,,y,), i = 
0, 1, ... ,s, also interpolates the couples (t,-,w;), i = l,...,r. That is, a(t) interpolates 
at * + 1 points, even though it has only degree s. If we define the abscissae 

{t <h < ■■■ <t k } = { Ci }U{Ti}, (1.13) 

and dispose them according to a Lobatto distribution in [0, 1] in order to get a formula 
of degree 2k, we have that 

1 k 

a'(t) T VH(a(t))dt = Y,bi<y , (t i ) T VH(a(t i )), (1.14) 

i=0 

and, consequently, the conservation condition becomes 

k 

Y,b i o'(t i ) T VH{o(t i ))=0, (1.15) 

;=0 



where, now, 



1 k t-t 



bi = / :— r*> 1 = 0,1,...,*. (1.16) 

j <>./-;'*' '/ 

The left-hand side of d 1 - 1 5b is called "discrete line integral" because, as will be 
clear in the sequel, the choice of the path a{t) is dictated by the numerical method 
by which we will solve problem dl.lt (see |[T3l for details). 

With these premises, in Section|H we devise such a method, able to fulfill ( 11.15b . 
after having set some preliminary results in Section [5] Before that, in Section |2] we 
state a few facts and notations concerning the shifted Legendre polynomials, which 
is the framework that we shall use to carry out the analysis of the methods. A few 
numerical tests are then reported in Section [5] and, finally, a few conclusions are 
given in Section [6] For sake of completeness, some properties of shifted Legendre 
polynomials are listed in the Appendix. 



6 



B. I. T. 



2 Preliminary results and notations 



The shifted Legendre polynomials, in the interval [0, 1], constitute a family of poly- 
nomials, {P„}„ g n> f° r which a number of known properties, named P1-P12, are re- 
ported in the Appendix. We now set some notations and results, to be used later. 
With reference to the abscissae ( 11.4b . let: 



\Pj(Cs) 



Pj 



j = 0,...,s, 



Pj e 



ps+l X j+l 



Jo l Pj(x)dx s 



I, 



> I; 



I, 



Remark 2 Observe that, from Pll, one obtains: 

h = 0. 

Furthermore, we set: 

V/ (I,, ...I, j •: R SXJ+l , Jj = (t ...tj) G M' s+lx -''+ 1 , 



/I 



\ 



n = 



and 



1 
1 

V 



o 



By virtue of P3 and P9, we deduce that 



J . . . J L» , 



and 



i = v&jGjDj 1 , i = x&PPi, J = 1,2, ... . 



(2.1) 

(2.2) 
(2.3) 

(2.4) 
(2.5) 
(2.6) 



(2.7) 



(2.8) 



(2.9) 
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Lemma 1 The matrix $ s — ( po 



/I 



T 

A-i o 



= ^G s = 



i+lxs+l - s nons i n g U i ar Moreover 



-1 

o '•• 



-1 

00 
10/ 



(2.10) 



a nonsingular matrix. 



with J s _\ € 1 

Proof & s is the transpose of the Gramian matrix defined by the linearly inde- 
pendent polynomials Po(c), . . . ,P s (c) at the distinct abscissae cq,...,c s and is, there- 
fore, nonsingular. The structure of J> s follows from (12.41 . The matrix J^_i is nonsin- 
gular since, from ( 12.101 i. & s is nonsingular, and rank(G. s ) = s. □ 



3 Matrix form of collocation methods 

In this section we deliberately do not care of the exactness of the discrete line integral, 
as stated by dl.141 . and in fact we choose k = s (and hence f; = c u i = 0, ...,s). We 
show that imposing the vanishing of the discrete line integral (condition ( 11.15b ) leads 
to the definition of the classical Lobatto IIIA methods. The reason why we consider 
this special situation is that the technique that we are going to exploit is easier to be 
explained, but at the same time is straightforwardly generalizable to the case k > s. 
As a by-product, we will gain more insight about the link between the new methods 
and the Lobatto IIIA class. For example, we will deduce that Lobatto IIIA methods 
(and, in general, all collocation methods) may be defined by means of a polynomial 
c(f) of degree not larger than that of the collocation polynomial u(t) (indeed, in the 
present case, degc(f) = degM(f) — 1). 

To begin with, let us consider the following expansion of d'(c): 

ff'(c) = £y^(c), (3.D 

where the (vector) coefficients jj are to be determined. Then, (I1.15l l becomes 

J-l s 

E WMQ)Vtf(c7(c,)) = 0, (3.2) 

7=0 i=0 

which will clearly hold true, provided that the following set of orthogonality condi- 
tions are satisfied: 

Yj = nj£biPj(ci)JVH(o(ct)), 7 = 0,. (3.3) 

1=0 
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where {rjj} are suitable scaling factors. We now impose that the polynomial 



s-l 



o-(c) =y + Y,Y) / Pj(x)dx 



satisfies (11.51 1. By setting 

/ 7o 
7= : 

one obtains (see (12.3l l-( |2~9"l i) 



ew, y 



( yl 



J s -\%h m f= ( ^ £? S G S D S 1 ) ®/2 m y=y-e0yo- 



Consequently, 



y= [2D s (,^ s G s y l {-el s )] ®h m y- 

On the other hand, the vector form of relations ( 13.31 ) reads 

y= (r^fl) ®/ 2m f, 

where T = diag(r/! , . . . , TJ S ) G R iXi and 

f (./„... ./.)' . ft = JVH(a(ci)), i = 0,...,s. 
Since F contains free parameters, we set 

r = Av. 



(3.4) 



(3.5) 



(3.6) 



(3.7) 



(3.8) 



(3.9) 



(3.10) 



Comparing (13.71 and ( 13.81 ). we arrive at the following block method, where now h 
denotes, in general, the used stepsize, 

A®I 2m y = hB®I Zm t, (3.11) 

with (see tilM ) 

A =(-«/,) , B = Q^C^LiflJ = (^i-iC^ifl) . (3.12) 

The following noticeable result holds true. 

Theorem 1 £ac/i row of the block method ( li./71 )-( li.72l ) defines a LMF of order s+ 1. 
77ie /flif row corresponds to the Lobatto quadrature formula of order 2s. 
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Proof For the first part of the proof, it suffices to show that the method is exact 
for polynomials of degree s + 1 . Clearly, it is exact for polynomials of degree 0, due 
to the form of the matrix A. We shall then prove that AJ? S = B£? s , that is (see (12. 3b , 



dOt , and dTT2l ). J s = B2> s . By virtue of (l3T2t . (1231 and (Ell, we have 

B# s = J S -\D S & T S _ X Q& S = J^_!D S [D7 1 0] = [S^ l s ) = J S1 



which completes the first part of the proof. For the second part, one has to show, by 
setting as usual e, the /th unit vector, that 



the vector containing the coefficients of the quadrature formula. From (13.12b , exploit- 
ing property P7 (see also ( 12.7b ). we obtain 



As an immediate consequence, the following result follows. 

Corollary 1 The block method A3.1 H - I[3.12]l collocates at the Lobatto abscissae ( 17.41 ) 
and has global order 2s. 

Proof The proof follows from known results about collocation methods (see, 
e.g., El Theorem II. 1.5]). □ 

Remark 3 In conclusion, the method corresponding to the pencil (A,B), as defined 
by ( 15.721 ), is nothing but the Lobatto IIIA method of order 2s. 

3.1 Link between a(c) and the collocation polynomial 

An important consequence of Theorem [T] and Corollary Q] is that the Lobatto IIIA 
method of order 2s may be also defined by means of an underlying polynomial, 
namely (7(c), of degree s instead of s + 1, as is the collocation polynomial associ- 
ated with the method ( 13.1 lb . 

The main aim of the present subsection is to elucidate the relation between these 
two polynomials. In what follows, we deliberately ignore the result obtained in The- 
orem[T]and CorollaryQ] so as to provide, among other things, an alternative proof of 
part of the statements they reportjj 

Let u(c) be the polynomial (11.6b (of degree 5+1) that collocates problem (11. lb at 
the abscissae (11.4b . The expansion of u'(c) along the shifted Legendre polynomials 
basis reads 



e^B=(b ... b s ), 



e T s B = lj3> l G l &l_ 1 Q 
= e[&U£l = (1 . 



1) 




(b ... b s ).D 



u 



'(c) = E CjPj(c). 



(3.13) 



3 The approach exploited in the proof of Theorem[T]turns out to be crucial to deduce the new methods 
presented in the next section. 
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Consequently, by setting 

/ go" 



gi=JVH(u(ci)), and C=(? 



( ( £° \ "\ 

V C / 



one obtains that ( 11.6b may be recast in matrix notation as ^ s <E> hm £ = g> or 

C = #7'«>^ m g. (3.14) 
We get the expression of m(c) by integrating both sides of (13.13b on the interval [0, c]: 



"(0 =yo+£ [ C Pj(x)dx + C s [ C p s (x)dx. (3.15) 

j = Q ''O 70 



By virtue of property Pll, we get 



s-1 



u{ci) =yo+Y / ^; i = 0, . . . , j. 

i~o ■><> 



(3.16) 



Setting z,- = k(c,-), i = 1, . . . ,s, z = (z\ , . . . ,z s ) T , and z = (jo,z r ) r , allows us to recast 
( 13.16b in matrix form. This is done by exploiting a similar argument used to get ( 13.61 ) 
starting from ( 13.4b : 



& s g s d: 



A®I 2m z = z-eOyo = ^-i®^2mC = 

^.g^d; 1 0]^®/ 2 „,t (3.17) 
Inserting (13. 14b into ( 13.171 ), and exploiting (12.81 ), yields 

= -^G s [D; 1 0] #7' ®/ 2 „,g = X -@ s G s & T s _ x £l ®b_ m g 

Thus, the collocation problem ( 11.6b defines the very same method arising from the 
polynomial a(c) (see (I3.llb - d3.12b ) with h=l. This implies that system ( 13.1 lb is 
a collocation method defined on the Lobatto abscissae c,, i = 0,. . . ,s, (therefore, a 
Lobatto IIIA method), and provides an alternative proof of Corollary Q] In particular, 
we deduce that 

u(ci)=yi = o(ci), / = 0,...,s. 
It follows that ( 13.15b becomes 



u{c) = a(c) + Cv / P s {x)dx 



(3.18) 
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and, after differentiating, 

u'(c) = a'{c) + ^P s (c). (3.19) 

We can obtain the expression of the unknown by imposing a collocation condition 
at any of the abscissae c, . For example, choosing c = c s = 1, yields 

Cv = u'{\) - a' (I) = f(y s ) - £ yj = fM - e T ®h,„y. (3.20) 

7=0 

This latter expression can be slightly simplified by considering that: 

(i) f(y s ) = (e T ®hn& which comes from the fact that the system SP^x = 

has solution x = e s + 1 (the nonsingularity of & s being assured by LemmaHJ; 

(ii) from (ETJl and d3~TTT i- (l3~T2l i, one has 

y = (D s , Q ) ® I 2 J = (D s &F_ , a & s 1 ) ® I 2 J 

= (d s {d; 1 Q)&- l )®i2j = (i s o)^" 1 ®i 2m i 

Thus, from (13 . 20t we get 

Cv = [{e T l)-(e T 0)] ®I 2 J = e T s+1 @-' ®/ 2m f. (3.21) 

The remaining collocation conditions, m'(c, ) = JVH(u(ci)),i = 0, ... ,s — 1, are clearly 
satisfied since the collocation polynomial u(c) is uniquely identified by the 5 + 2 lin- 
early independent conditions in ( 11.61 ). Nonetheless, they can be easily checked after 
observing that, from ( 13.18b . (ii), and ( 13.21b , 

1= ( J) =0>J x ®1tJ. 



Therefore, from ( I3.13l l, ( 13.14b . and (13.18b . one obtains, 
(u'(c )\ 

u' = : = & s ®7 2m C = ® 7 2m f = f. 

\ "'(*)/ 

That is (see d3~9V), m'(q) = JVH(u( Ci )), i = Q,...,s. 



4 Derivation of the methods 

The arguments for deriving the methods in the general case where k > s (which makes 
the discrete line integral exact) are a straightforward extension of what stated above. 
In particular, let us consider again the expansion ( 13.1b - (13.4b of the polynomial <7(c). 
Then, condition ( 11.15b can be recast as (compare with (13.2b ) 



s— 1 k 

£yJ£fc// J / -(^)VH(c7(? i )) = o, 



(4.1) 
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which will clearly hold true, provided that the following set of orthogonality condi- 
tions are satisfied (compare with (13.31 l. see also (11.161 ): 



Yj = flj£,b i P j (t i )JVH(G(t i )), j = 0,...,s-l, 



(4.2) 



;=o 



where are suitable scaling factors. According to ( 13.10b . we choose them as 
T]j = 2j + 1, j = 0, . . . , s — 1 . The vector y (see ( 13.51 )) is then obtained by imposing that 
the polynomial c(c) in (13. 4b satisfies the interpolation constrains ( 11. 5b and ( 11.121 ). 
In so doing, one obtains a block method characterized by the pencil (A,B), where 
the two k x k + 1 matrices A and B are defined as follows. In order to simplify the 
notation, we shall use a "Matlab-like" notation: let ind s £ W + and ind r £ W be the 
vectors whose entries are the indexes of the main abscissae cq < ■ ■ ■ < c s in ( 11.4b and 
of the silent ones Ti < • • • < T r in dl.l lb . respectively, within the Lobatto abscissae 
fo < • • • < /fc, as defined in dl . 1 3b . Then, the orthogonality conditions (14.2b will define 
the first s rows of A and (compare with (13.12b ): 

A(l : s,ind s ) = ( -e I s ) , B{\ :*,:) = (^D S ^ T £2 

where (see (t23l— d2T6b and ( fTUb ) 

M('o) •■ 



(4.3) 



and (see ( 11.7b ) 



Q 



\Po(tk) 
\ 




pir+l Xi 



(4.4) 



fk+lxk+l 



(4.5) 



bk 



On the other hand, the interpolation conditions for the silent stages ( 11.12b define the 
last r rows of the matrix A (the corresponding rows of B are obviously zero): 



A(s + 1 : k, ind r ) = Ir, 

A(s+ 1 : k, ind s ) = -J r ( -e I„)] - , 

where I r is the identity matrix of dimension r, e=(l,...,l) r 
vector (of dimension 5+1), and 

/./,; /Uv;i/.v... ^P s _ x {x)dx 



(4.6) 



I ■ 



e\ is the first unit 



\J Q Tr P (x)dx ... .^P^^x)^, 



The following result generalizes Theorem[T]to the present setting (the proof being 
similar). 



4 As a further convention, the entries not explicitly set are assumed to be 0. 
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Theorem 2 Each row of the block method A4.3[ — fcL6i defines a LMF of order at least 
s. The s-th row corresponds to the Lobatto quadrature formula of order 2k. 

Definition 1 We call the method defined by the pencil (A,B) i njA3\ - l4.6\ a "Hamil- 
tonian B VM with k steps and degree s", hereafter HBVM(£,s)|j 

Remark 4 The structure of the nonlinear system associated with the HBVM(k,s) is 
better visualized by performing a permutation of the stages that splits, into two block 
sub-vectors, the fundamental stages and the silent ones. More precisely, the permuted 
vector of stages, say z, is required to be: 



= bo,yi r--J.,;^l ,W2y 
fundamental stages silent stages 



This is accomplished by introducing the permutation matrices W £ R and W\ G 
R *+lx*+l > suchthat 



( 2 



W 



\k+l 
It is easy to realize that 



ind s (2:s+l) 
ind r 




ind v 



WAW, 



-e h 0. sxr 
-ao -A\ I, 



WBW, 



b B, B 2 

O r xs Orxr 



where [—ao, — A\\ coincides with A(s+ 1 : k,md s ) in i4.6i , while [bo, B\ , B2] matches 
the matrix B(l : s,:) in (14.3b . The HBVM(k,s) then takes the form: 



-e I s (W 
-a -A\ I r 



®h m z = h 



b B } B 2 

O r xs O r xr 



)JVH(z). 



(4.7) 



The presence of the null blocks in the lower part ofWBWf clearly suggests that the 
(generally nonlinear) system ( 14.7b of (block) size k is actually equivalent to a system 
having (block) size s. Indeed, we can easily remove the silent stages, 

w = a ®yo+A\ ®h m y, 

and obtain 



y = e®yo + hbo®{Tm(yo)) + hBi®JVH{y) 

+hB 2 ®JVH{ao®yo+Ai®h m y). 



(4.8) 



We refer to |2f for an alternative technique to reduce the dimension of system d4.7K 



Indeed, the pencil (A,B) perfectly fits the framework of block BVMs (see, e.g., (3)). 
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Remark 5 As was shown in the previous section, when k = s, the HBVM(s,s) coin- 
cides with the Lobatto III A method of order 2s. More in general, fork>s, by summing 
up ( I4.3D - R61 ), we can cast HBVM(k,s) as a Runge-Kutta method with the following 
tableau: 



to 






JD S 3P T Q. 


k 






h) ... b k 



(4.9) 



where 

/J ,0 P Q (x)dx ... $P s ^(x)dx\ 

y=\ \ \ et i+lxs . 

\$P (x)dx... J^P s ^(x)dx) 
We observe that the k+l x k + 1 matrix 

C = JD S & T Q. (4.10) 

appearing in ( 14.91 ) has rank s, thus confirming that the computational cost per itera- 
tion depends on s, rather than on k (see fiZS for more details and a practical example 
of Butcher tableau concerning the method HBVM(6,2)). 

By the way, we observe that, when s = 1, HBVM(k, 1 ) are nothing but the "s-stage 
trapezoidal methods ", defined in HI Oil , based on the Lobatto abcissae. In such a case, 
the matrix C becomes 

C") 

' (f>o ... h). 

Similarly, for s = 2 and k = 4, HBVM(4,2) coincides with the fourth-order method 
presented in M3\ Section 4.2], able to preserve polynomial Hamiltonians of degree 
four. 

Concerning the order of convergence, the following result generalizes that of 
Corollary [1] 

Corollary 2 The HBVM(k,s) ( !4.3l )-( l4T6t has order of convergence p = 2s. 

Proof By virtue of Theorem |2 the corresponding Runge-Kutta method ( 14.91 ) sat- 
isfies the usual simplifying assumptions B(2k) and C(s). If we are able to prove 
D(s — 1), from the classical result of Buthcher (see, e.g., J9) Theorem 5.1]), it will fol- 
low that the method has order p = 2s. With reference to J4.91 , the condition D(s — 1) 
can be cast in matrix form, by introducing the vectors e = (1,... ,l) r G W~ l , e = 
(1, . . . , l) r S R* +1 , and the matrices 

G = diag(l,...,5-1), r = diag(fo,....fc), V = (/jt 1 1 )GR* +lx '- 1 , 



as 



QV T Q {JD S & T Q) = (ee T -V T T)Q, 
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i.e., 



S?D S J T Q.VQ = (ee T - TV) 



(4.11) 



Since the quadrature is exact for polynomials of degree 2s — 1 < 2k— 1, one has 



(J^QVQ) = ( [" ' Pi-i{x)dx{jti 



,e=o 



1 rt 



JO 



P,-i{x)Ax{jt'- l )At 



S n -y^ P i - l (x)x J dxj , i=l,...,s, j = l,...,s-l, 

where the last equality is obtained by integrating by parts, with <5,i the Kronecker 
symbol. Consequently, 

{&D s J T QVQ) ij = (l-ZjltPth) £Pi{x)x j dx^j = (1 -</_!), 

i=l,...,k+l, j = l,...,s-l, 
that is, (14.11b . where the last equality follows from the fact that 

V TfeJ%(f) / P e (x) X jdx =tJ, j = 1, ... ,S- 1. □ 

/= ■/<> 



An additional, remarkable property of such methods is gained, provided that the 
abscissae {to, ■ ■ ■ ,?<.} ( 11.131 ) are symmetrically distributed (as is the case of the Lo- 
batto abscissae here considered). For this purpose, we need to introduce some nota- 
tions and preliminary results. Let us define the matrix 



/ 1\ 



E n = 



which, when applied to a vector of length n, reverses the order of its entries. We also 
set 



1 1 



1 1 



/ 

re 

V 



\ 



(-irv 



(4.12) 



The following preliminary result holds true. 

Lemma 2 If the abscissae M.13\ are symmetric, then matrix H4.10\l satisfies: 
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Proof From the symmetry of the abscissae it easily follows that (see d 1 - 1 6b and 

E k+ iQE k+l = £2. 
From property P6, we have that (see (|4.4| i) 

& T E k+x = F& T . 

Moreover, by considering that (see fll.4| >) 

/ %P (x)dx ... %P s -i(x)dx \ 
LJ= : : 

\Jl i P (x)dx... f!; I\ { {x)dx) 

again from P6 we see that 

E S LJ =LJF. 

Finally, from ( 14.101 ) we obtain 

E k LCE k+ \ = 

= (E k LJ?)D s {^ T E k+l ){E k+l QE k+1 ) 

= LJ ; FD S F£? T Q.=LJ?D S Z? >T Q.=LC. □ 

As a consequence, we have the following result. 

Theorem 3 If the abscissae ( 17.731 ) are symmetric, then the method H4.3]) - fc4~6]) (i.e., 
H4.9\l ) is symmetric, that is, it is self-adjoint. 

Proof Indeed, the discrete solution, y, satisfies the equation (see ( |4.9b — (14. 10b 
and gJ3) 

L®hmf = hLC®I lm f($). 

Considering that E k LE k+ \ = — L and, from Lemma [2] E k LCE k+ \ = LC, one then 
obtains 

L®hm{Ek+i®hmt) = -hLC®I 2m {E k+ i®hmf{y)) 
= -hLC®hmf(E k+l ®hmy)- 

The thesis then follows by observing that the vector E k+ \ ®h m y contains the time- 
reversed discrete solution. □ 

The next theorem summarizes the results about HBVM(A:,.s). 

Theorem 4 (Main Result) Foralls = 1,2,..., andk>s, the HBVM(k,s) method: 

1. has order of accuracy 2s; 

2. is energy-preserving for polynomial Hamiltonians of degree not larger than 2k / s; 

3. for general C^ 2k+l ' 1 Hamiltonians, the energy error at each integration step is 
0(h 2k+l ), ifh is the used stepsize;^ 

6 Consequently, on any finite interval the global energy error is not larger than 0(h 2k ). 
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4. is symmetric and, therefore, precisely A-stable. 

Proof Item 1 follows from Corollary [2] Item 2 follows from the fact that, for 
such polynomial Hamiltonians, the vanishing discrete line integral equals the con- 
tinuous line integral (see ( 1 1.141 ) and ( 11.15b ). Similarly, Item 3 follows from the fact 
that, by using arguments similar to those used in Remark Q] (see (11.81 )). the energy 
error per integration step equals the quadrature error of the Gauss-Lobatto formula 
of order 2k. Finally, Item 4 follows from Theorem[3] since the Lobatto abscissae {f, } 
are symmetrically distributed. □ 

Remark 6 From the result of Theorem^ we can then concude that HBVM(k,s) is 
optimal, both from the point of view of the order and stability properties. Moreover, 
its computational cost, as observed in Re?narks^\and\5\ is seen to depend on s, rather 
than on k. 

5 Numerical Tests 

We here report a few numerical tests, in order to show the potentialities of HB VM(k, s) . 

Let then consider, at first, the Hamiltonian problem characterized by the polyno- 
mial Hamiltonian (4.1) in Q, 

ff (^) = y-2 + 30 + T-y + 6' (5 - 1} 

having degree v = 6, starting at the initial point yo = (q(Q), p(0)) T = (0, 1) T . For such 
a problem, in [7| it has been experienced a numerical drift in the discrete Hamiltonian, 
when using the fourth-order Lobatto IIIA method^ with stepsize h = 0.16. This is 
confirmed by the plot in Figure lBTTl where a linear drift in the numerical Hamiltonian 
is clearly observable. On the other hand, by using the fourth-order HBVM(6,2) with 
the same stepsize, the drift disappears, as shown in Figure 15.21 since such method 
exactly preserves polynomial Hamiltonians of degree up to 6. Moreover, the order 
of convergence p = 4 is (numerically) confirmed by the results listed in Table 15.11 
where the used stepsizes h, the maximum estimated error (obtained as the difference 
of two consecutive solutions), and the estimated order of convergence are listed. 

The second test problem is the Fermi-Pasta-Ulam problem (see JU Section 1.5.1]), 
defined by the Hamiltonian 

j m m m 

H(p,q) = - £ + Pli) + -r E (?2i - f + £ [qu+\ - quf , (5.2) 

1 !=1 4 1=1 1=0 

with qa = qom+\ = 0, m = 3, (O = 50, and starting vector 

Pi = 0, <?, = (/-l)/10, i=l,...,6. 

In such a case, the Hamiltonian function is a polynomial of degree 4, so that the 
fourth-order HBVM(4,2) method, which is used with stepsize h = 0.05, is able to 



7 Such method coincides with the HBVM(2,2) above described. 




Fig. S.2 Fourth-order HBVM(6,2) method, h = 0. 16, problem l5l\ . 
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Table 5.1 Numerical order of convergence for the HBVM(6,2) method, problem <5~Tl 



/; 


0.32 0.16 0.08 0.04 0.02 


error 


2.288 10" 2 1.487- 10" 3 9.398 • 10 -3 5.890 - lO"" 3.684- 10" ' 


order 


3.94 3.98 4.00 4.00 


Table 5.2 Numerical order of convergence for the HBVM(4,2) method, problem 15.2). 


/; 


1.6- 10"^ 8-10" J 4-10~ J 2-10- J Kr j 




error 


3.030 1.967- Kr 1 1.240- 10~ 2 7.761 -10 -4 4.853- lO -3 


order 


3.97 3.99 4.00 4.00 


Table 5.3 Numerical order of convergence for the HBVM(6,2) method, problem 15.3). 


/; 


3.2- 10~ 2 1.6- 10~ 2 8-10 _J 4-Kr J 210" J 


error 


3.944 -lO" 6 2.635 -lO" 7 1.729 lO" 8 1.094- 10- tJ 6.838- 10" 11 


order 


3.90 3.93 3.98 4.00 



exactly preserve the Hamiltonian, as confirmed by the plot in Figure 15741 whereas the 
fourth-order Lobatto IIIA method provides the result plotted in Figure fOl Moreover, 
in Table 15.21 we list corresponding results as in Table l5.ll again confirming the fourth- 
order convergence. 

In the previous examples, the Hamiltonian function was a polynomial. Neverthe- 
less, as is easily argued from Theorem[4] HBVM(fc,s) are expected to produce a prac- 
tical conservation of the energy when applied to systems defined by a non-polynomial 
Hamiltonian function which are sufficiently differentiable. As an example, we con- 
sider the motion of a charged particle in a magnetic field with Biot-Savart potential^ 
It is defined by the Hamiltonian 



H(x,y,z,x,y,z) = (5.3) 
1 



2m 



. 2 / \ 2 



with p = y/ 'x 2 +y 2 , (X = eBo, m is the particle mass, e is its charge, and Bo is the 
magnetic field intensity. We have used the values 

m = 1 , e = — 1, Bq = 1 , 

with starting point 

x = 0.5, y = \0. z = 0, ±=-0.1, y = -0.3, z = Q. 

In Figure 1531 the trajectory of the particle in the interval [0, 10 3 ] is plotted in the 
phase space. As one can see, it is a helix that wings downward. By using the fourth- 
order Lobatto IIIA method with stepsize h = 0. 1 , a drift in the numerical Hamiltonian 
can be again observed (see Figure IBToT i. so that the method does introduce a friction. 
When using the HBVM(4,2) method with the same stepsize, the drift disappears and 
the Hamiltonian turns out to be almost preserved (see Figure [57Tb . As expected, the 



As an example, this kind of motion causes the well known phenomenon of aurora borealis. 
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Fig. 5.4 Fourth-order HBVM(4,2) method, h = 0.05, problem <BT2l 
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result improves if we increase k: the plot in Figure l5~8l has been obtained by using the 
HBVM(6,2), from which one realizes that a practical preservation of the Hamiltonian 
is reached. Finally, the data listed in Table 15 . 31 confirm the fourth-order convergence 
of the latter method. 



6 Conclusions 

In this paper a new class of numerical methods, able to preserve polynomial Hamil- 
tonians, has been studied in details. From the analysis, it turns out that such methods 
can be regarded as a generalization of collocation Runge-Kutta Lobatto IIIA meth- 
ods. Nevertheless, the fact of being characterized by a matrix pencil, perfectly fits 
the framework of block BVMs, so that we have named them Hamiltonian BVMs 
(HBVMs). A number of numerical tests prove their effectiveness in preserving the 
Hamiltonian function when evaluated along the numerical solution, as well as con- 
firm the predicted order of convergence. Possible different choices of the abscissae, 
as well as the actual efficient implementation of the methods, will be the subject of 
future investigations. 
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Appendix: some properties of shifted Legendre polynomials 

A number of useful properties of shifted Legendre polynomials are here summarized: 
for their proof see any book on special functions (e.g., [fl~)). 

PI. Generalized Rodrigues formula: for all n = 0, 1 , . . . , P„(x) has degree n and can 
be defined as 

P2. Lobatto quadrature: the Lobatto abscissae {c,} dl.4t , of the formula of degree 
2s, are the zeros of the polynomial 

(x 2 -4P s '(x), 

where P^(x) denotes the derivative of P s (x). The corresponding weights (11.7b are 
given by: 

which are, therefore, all positive. 




Fig. 5.6 Fourth-order Lobatto IIIA method, h = 0.1, problem \53\ . 
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P3. Orthogonality: 

f 1 1 

/ P n (x)P m (x)dx= - — — 8 nm , n = 0, 1 , . . . , 
Jo 2n + 1 

where, as usual, 8„„, denotes the Kronecker delta. 
P4. Recurrence formula: by setting hereafter P-\{x) = and Po{x) = 1, 

(n + 1 )P n+l (x) = (2n + l)(2x - 1 )P n {x) - nP„^(x) , « = 0, 1 , . . . . 

P5. Explicit formula: 

ft(x) = (-l)-t(")( n 1"'') (-*)', n = 0,l,.... 
P6. Symmetry: 

P„(l-x) = (-l)"P„(x), « = 0,1,.... 
P7. Symmetry at the end-points: 

P„(0) = (-l) n , P„(l) = l, « = 0,1,.... 

P8. Derivatives: 

2(2n + 1 )P„(x) = — [P n+i (x) - P n - 1 (*)] , « = 0, 1 , . . . . 

P9. Integrals: 

2/ Po(t)dt = 2x = P l (x)+P (x), 



2(2« + l)/ P n {t)dt =P n+l (x)-P n - 1 (x), n = l,2,.... 
Jo 

P10. Shifted Legendre differential equations. The shifted Legendre polynomials sat- 
isfy the second order differential equation: 

— [(x 2 - x)P' n {x)} +n(n+ l)P n {x)=0, « = 0,1,.... 



Pll. From P2 and P10, it follows that, if dl.4t are the Lobatto abscissae of the for- 
mula of order 2s (i.e., exact for polynomials of degree 2s — 1), then 



P12. A few examples: 



P s (x)dx = 0, i = 0,1,..., s. 



P (x) EE 1, 

Pi{x) =2x-l, 

P 2 (x) = 6x 2 -6x+l, 

P 3 (x) = 20x 3 - 30x 2 + I2x - 1 , 
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